function [C] = binopalonglastdim(A, B, hBinOps, doSqueeze)
%BINOPALONGLASTDIM 沿最后一维进行二元运算
%
% Input Arguments
% # A: 数组
% # B: 数组，最后一维尺寸与A最后一维尺寸相同
% # hBinOps: 二元操作符句柄，默认为mtimes
% # doSqueeze: 逻辑标量，true表示移除结果中尺寸为1的维度，默认为true
%
% Output Arguments
% # C: 数组，最后一维尺寸与A最后一维尺寸相同，除最后一维外的尺寸由hBinOps结果决定
%
% Assumptions and Limitations
% # 假设MATLAB使用列优先索引

% NOTE: 本函数使用线性索引
if nargin < 3
    hBinOps = @mtimes;
end
if nargin < 4
    doSqueeze = true;
end

szA = size(A);
szB = size(B);
assert(szA(end) == szB(end));
C1 = hBinOps(lastdimsubarray(A, 1), lastdimsubarray(B, 1));
szC1 = size(C1);
nElemC1 = prod(szC1);
C = NaN([szC1 szA(end)]);
C(1:nElemC1) = C1;
for i=2:szA(end)
    linearIndC = ((i-1)*nElemC1+1):(i*nElemC1);
    C(linearIndC) = hBinOps(lastdimsubarray(A, i), lastdimsubarray(B, i));
end
if doSqueeze
    C = squeeze(C);
end
end